       m_a   = (*fields_ptr).u[M_A  ]GRD_IND;
       m_a_t = (*fields_ptr).u[M_A_T]GRD_IND;
       m_a_x = (*fields_ptr).u[M_A_X]GRD_IND;
       m_a_y = (*fields_ptr).u[M_A_Y]GRD_IND;

       m_b   = (*fields_ptr).u[M_B  ]GRD_IND;
       m_b_t = (*fields_ptr).u[M_B_T]GRD_IND;
       m_b_x = (*fields_ptr).u[M_B_X]GRD_IND;
       m_b_y = (*fields_ptr).u[M_B_Y]GRD_IND;

       m_c   = (*fields_ptr).u[M_C  ]GRD_IND;
       m_c_t = (*fields_ptr).u[M_C_T]GRD_IND;
       m_c_x = (*fields_ptr).u[M_C_X]GRD_IND;
       m_c_y = (*fields_ptr).u[M_C_Y]GRD_IND;

	norm_m = sqrt(m_a*m_a + m_b*m_b + m_c*m_c);
        m_a = m_a/norm_m;
        m_b = m_b/norm_m;
        m_c = m_c/norm_m;


		g_up_xx = auxfields.u_aux[G_UP_XX]GRD_IND;
		g_up_xy = auxfields.u_aux[G_UP_XY]GRD_IND;
       	g_up_yy = auxfields.u_aux[G_UP_YY]GRD_IND;

//       	Sdet = sqrt(g_up_xx*g_up_yy - g_up_xy*g_up_xy);
       
        n_up_x = g_up_xx*nx + g_up_xy*ny;
        n_up_y = g_up_xy*nx + g_up_yy*ny;

        norm_n = sqrt(n_up_x*nx + n_up_y*ny);

// Note we don't normalize the normals! At the end you see that the normalization drops out (because the eigenvalue is proportional to the normal), so is really irrelevant.


 SP = norm_n;

nV_a = n_up_x*(m_a_x - m_a_x_n) + n_up_y*(m_a_y - m_a_y_n);
nV_b = n_up_x*(m_b_x - m_b_x_n) + n_up_y*(m_b_y - m_b_y_n);
nV_c = n_up_x*(m_c_x - m_c_x_n) + n_up_y*(m_c_y - m_c_y_n);

V_a = 0.5*((m_a_t - m_a_t_n) + nV_a/norm_n)/norm_n; 
V_b = 0.5*((m_b_t - m_b_t_n) + nV_b/norm_n)/norm_n; 
V_c = 0.5*((m_c_t - m_c_t_n) + nV_c/norm_n)/norm_n; 


//#ifdef FULL_PROJECTION

//V_dot_m = m_a*V_a + m_b*V_b + m_c*V_c;

//V_a = V_a - m_a*V_dot_m;
//V_b = V_b - m_b*V_dot_m;
//V_c = V_c - m_c*V_dot_m;

//#endif


	     (*derivs_ptr).u[M_A_T]GRD_IND = (*derivs_ptr).u[M_A_T]GRD_IND 
			                   - Tau_p*SP*norm_n*V_a;
	     (*derivs_ptr).u[M_A_X]GRD_IND = (*derivs_ptr).u[M_A_X]GRD_IND 
					   - Tau_p*SP*nx*V_a;
	     (*derivs_ptr).u[M_A_Y]GRD_IND = (*derivs_ptr).u[M_A_Y]GRD_IND 
					   - Tau_p*SP*ny*V_a;

	     (*derivs_ptr).u[M_B_T]GRD_IND = (*derivs_ptr).u[M_B_T]GRD_IND 
			                   - Tau_p*SP*norm_n*V_b;
	     (*derivs_ptr).u[M_B_X]GRD_IND = (*derivs_ptr).u[M_B_X]GRD_IND 
					   - Tau_p*SP*nx*V_b;
	     (*derivs_ptr).u[M_B_Y]GRD_IND = (*derivs_ptr).u[M_B_Y]GRD_IND 
					   - Tau_p*SP*ny*V_b;

	     (*derivs_ptr).u[M_C_T]GRD_IND = (*derivs_ptr).u[M_C_T]GRD_IND 
			                   - Tau_p*SP*norm_n*V_c;
	     (*derivs_ptr).u[M_C_X]GRD_IND = (*derivs_ptr).u[M_C_X]GRD_IND 
					   - Tau_p*SP*nx*V_c;
	     (*derivs_ptr).u[M_C_Y]GRD_IND = (*derivs_ptr).u[M_C_Y]GRD_IND 
					   - Tau_p*SP*ny*V_c;
